require(colorout)
require(sqldf)
require(arm)
require(maps)
require(mapproj)
require(maptools)
require(sp)
require(rgdal)
require(RPostgreSQL)

remove(list=objects())
options(digits=2, scipen=9, width=110, java.parameters = "-Xrs")

################################################################################################################
# helper functions

username <- system("echo $USER", intern=TRUE)
conn <- dbConnect(dbDriver("PostgreSQL"), user=username, password="", host="localhost", port=5432, dbname=username)
SQL <- function(x) 
  dbGetQuery(conn, x)
SQL_EXECUTE <- function(x)
  system(paste0("psql -d $USER -U $USER -c '", x, "';"))

ColConvert <- function(x, rng=c(-0.2, 0.2)) {
  x <- pmax(rng[1], pmin(rng[2], x))
  mu <- mean(rng)
  x <- x - mu
  rng <- rng - mu
  rng <- abs(rng)^0.5 * sign(rng)
  x <- abs(x)^0.5 * sign(x)
  x <- round(1000 * (x - rng[1]) / (rng[2] - rng[1])) + 1
  return(x)
}

CMUSansSerif <- Type1Font(
  "CMUSansSerif", 
  c("fonts/cmuss1.afm",
    "fonts/cmuss2.afm",
    "fonts/cmuss3.afm",
    "fonts/cmuss4.afm"))
postscriptFonts(CMUSansSerif=CMUSansSerif)
pdfFonts(CMUSansSerif=CMUSansSerif)

################################################################################################################
# get data from database

dat <- SQL(paste0("
  select 
  fips, 
  case when race = 'Other' then 'Asian/Other' else race end as race, 
  gender, 
  count(*) as n, 
  avg(yhat) as yhat
  from poststrat a 
  inner join final b on (a.randomid=b.randomid)
  where voted2012 = 1
  group by 1, 2, 3
  order by 1, 2, 3;
"))
options(sqldf.driver = "SQLite")

################################################################################################################
# load shapefiles

shp <- list()
shp$state <- readShapePoly(
  "data/st99_d00_shp/st99_d00.shp", 
   proj4string=CRS("+proj=longlat"))
shp$county <- readShapePoly(
  "data/co99_d00_shp/co99_d00.shp", 
   proj4string=CRS("+proj=longlat"))

projection <- "+proj=lcc +lat_0=39 +lon_0=-96"
for (i in 1:length(shp))
  shp[[i]] <- spTransform(shp[[i]], CRS(projection))

# county centroids, approximate from shapefile polygons
shp$county@data$fips <- paste0(shp$county@data$STATE, shp$county@data$COUNTY)
tmp <- getSpPPolygonsLabptSlots(shp$county)
shp$county@data$x <- tmp[,1]
shp$county@data$y <- tmp[,2]
tmp <- shp$county@data
centroids <- sqldf("select fips, avg(x) as x, avg(y) as y from tmp group by 1")

################################################################################################################
# figure 5

agg <- sqldf(paste0("
  select 
  fips, race, 
  sum(n) as n, 
  sum(yhat * n) / sum(n) as yhat
  from dat
  group by 1, 2
"))
cols <- paste0(colorRampPalette(c("dark red", "grey85", "dark blue"))(1001), "BF")
agg$col <- cols[ColConvert(agg$yhat, rng=c(0.25, 0.75))]

widths <- c(5,5,1.2)
heights <- c(3.5,3.5)
pdf("figures/1percent_figure5.pdf", width=9, height=9 * sum(heights) / sum(widths), family=CMUSansSerif)
par(cex.main=1.25)

nf <- layout(
  matrix(c(01, 02, 05,
           03, 04, 05), nrow=2, byrow=TRUE),
  widths=widths,
  heights=heights,
  respect=TRUE)
#layout.show(nf)

par(mar=c(0,0,0,0))
for (race in c("White", "Black", "Hispanic", "Asian/Other")) {
  tmp <- merge(x=agg[agg$race == race,], y=centroids)
  plot(shp$state, xlim=c(-2.25e6, 2.15e6), ylim=c(-0.9e6, 0.8e6), border="dark grey", lwd=0.5)
  symbols(tmp$x, tmp$y, 
    circles=3000 * sqrt(tmp$n), inch=FALSE, 
    add=TRUE, fg=NA, bg=tmp$col)
  title(paste0("\n\n\n", race))
}

par(mar=c(0,0,0,0), lty=0)
plot(0, 0, type="n", xlab="", ylab="", xlim=c(0,1), ylim=c(0,1), axes=F)
plotrix::color.legend(
  0.15, 0.025, 0.35, 0.975, 
  rect.col=substr(cols[ColConvert(seq(0.25,0.75,length.out=101), rng=c(0.25,0.75))], 1, 7), 
  gradient="y", align="rb", 
  legend=c("  25%\n  Obama", "  50%\n  Obama", "  75%\n  Obama"),
  cex=0.8)
par(lty=1)

dev.off()

################################################################################################################
# figure 6

agg <- sqldf(paste0("
  select 
  fips, 
  'White ' || gender as gender, 
  sum(n) as n, 
  sum(yhat * n) / sum(n) as yhat
  from dat
  where race = 'White' and gender in ('Female', 'Male')
  group by 1, 2
"))
cols <- paste0(colorRampPalette(c("dark red", "grey85", "dark blue"))(1001), "BF")
agg$col <- cols[ColConvert(agg$y, rng=c(0.25, 0.75))]

gap <- sqldf(paste0("
  select 
  fips, 
  sum(n) as n, 
  avg(case when gender = 'Female' then yhat else NULL end) as female, 
  avg(case when gender = 'Male' then yhat else NULL end) as male
  from dat
  where race = 'White' and gender in ('Female', 'Male')
  group by 1
"))
gap$gap <- pmax(0, gap$female - gap$male)
cols <- paste0(colorRampPalette(c(rgb(0.5, 0, 0.5), "grey85", rgb(0, 0.5, 0)))(1001), "BF")
gap$col <- cols[ColConvert(gap$gap, rng=c(0, 0.15))]
gap <- merge(x=gap, y=centroids)

widths <- c(5,1)
heights <- rep(3.5,6)/2
pdf("figures/1percent_figure6.pdf", width=5.5, height=5.5*sum(heights)/sum(widths), family=CMUSansSerif)
par(cex.main=1.25)

nf <- layout(
  matrix(c(01, 04,
           01, 04,
           02, 04,
           02, 05,
           03, 05,
           03, 05), ncol=2, byrow=TRUE),
  widths=widths,
  heights=heights,
  respect=TRUE)
#layout.show(nf)

par(mar=c(0,0,0,0))
for (gender in c("White Female", "White Male")) {
  tmp <- merge(x=agg[agg$gender == gender,], y=centroids)
  plot(shp$state, xlim=c(-2.25e6, 2.15e6), ylim=c(-0.9e6, 0.8e6), border="dark grey", lwd=0.5)
  symbols(tmp$x, tmp$y, 
    circles=3000 * sqrt(tmp$n), inch=FALSE, 
    add=TRUE, fg=NA, bg=tmp$col)
  title(paste0("\n\n\n", gender))
}

plot(shp$state, xlim=c(-2.25e6, 2.15e6), ylim=c(-0.9e6, 0.8e6), border="dark grey", lwd=0.5)
symbols(gap$x, gap$y, 
  circles=3000 * sqrt(gap$n), inch=FALSE, 
  add=TRUE, fg=NA, bg=gap$col)
title(paste0("\n\n\n", "White Gender Gap"))

par(mar=c(0,0,0,0), lty=0)
plot(0, 0, type="n", xlab="", ylab="", xlim=c(0,1), ylim=c(0,1), axes=F)
cols <- paste0(colorRampPalette(c("dark red", "grey85", "dark blue"))(1001), "BF")
plotrix::color.legend(
  0.15, 0.025, 0.35, 0.975, 
  rect.col=substr(cols[ColConvert(seq(0.25,0.75,length.out=101), rng=c(0.25,0.75))], 1, 7), 
  gradient="y", align="rb", 
  legend=c("  25%\n  Obama", "  50%\n  Obama", "  75%\n  Obama"),
  cex=0.8)

plot(0, 0, type="n", xlab="", ylab="", xlim=c(0,1), ylim=c(0,1), axes=F)
cols <- paste0(colorRampPalette(c(rgb(0.5,0,0.5), "grey85", rgb(0,0.5,0)))(1001), "BF")
plotrix::color.legend(
  0.15, 0.025, 0.35, 0.975, 
  rect.col=substr(cols[ColConvert(seq(0.25,0.75,length.out=101), rng=c(0.25,0.75))], 1, 7), 
  gradient="y", align="rb", 
  legend=rev(c("  Female\n  +15%\n  Obama", "  Female\n  +7.5%\n  Obama", "  Female\n  +0%\n  Obama")),
  cex=0.8)

dev.off()

#########################################################################################################
# figure 7

pdf("figures/1percent_figure7.pdf", width=6.5, height=3.75, family=CMUSansSerif)
par(mfrow=c(1,2), mar=c(2,2,2.5,1), mgp=c(1,0.01,0), tck=0.005, lwd=0.5, pty="s", 
  cex.main=0.85, cex.lab=0.75, cex.axis=0.65)

agg <- sqldf("
  select fips, 
  sum(n) as n, 
  sum(case when race = 'Hispanic' then n else 0 end) / sum(n) as x, 
  sum(case when race = 'Hispanic' then n * yhat else 0 end) / sum(case when race = 'Hispanic' then n else 0 end) as y
  from dat
  group by 1
")
agg <- agg[complete.cases(agg),]
agg$x <- log10(agg$x)

plot(0, 0, type="n", axes=FALSE,
  xlab="Percent Hispanic (Log Scale)", ylab="Hispanic Obama Support", 
  ylim=c(0,1), xlim=c(log10(0.0001), log10(1)),
  main="Higher Hispanic Obama Support\nin Majority Hispanic Areas")
symbols(agg$x, agg$y, circles=sqrt(agg$n) / 500, 
  bg=rgb(0.5, 0.5, 0.5, alpha=0.1), fg=rgb(0.5, 0.5, 0.5, alpha=0.75), 
  add=TRUE, inch=FALSE)
L <- loess(y ~ x, weights=n, span=1, data=agg)
ord <- order(agg$x)
points(agg$x[ord], L$fit[ord], type="l")
axis(1, at=-4:0, lab=paste0(formatC(100 * 10^(-4:0)), "%"), lwd=par()$lwd)
axis(2, at=seq(0,1,by=0.25), lab=paste0(100*seq(0,1,by=0.25), "%"), lwd=par()$lwd)
box(lwd=par()$lwd)

######

agg <- sqldf("
  select 
  fips, 
  sum(n) as n, 
  sum(n * yhat) / sum(n) as x, 
  sum(case when gender = 'Female' then yhat * n else 0 end) / sum(case when gender = 'Female' then n else 0 end) - 
    sum(case when gender = 'Male' then yhat * n else 0 end) / sum(case when gender = 'Male' then n else 0 end) as y 
  from dat
  where race = 'White'
  group by 1")
agg$y <- pmin(0.15, pmax(0, agg$y))
agg <- agg[complete.cases(agg),]

plot(0, 0, type="n", axes=F,
     ylab="White Gender Gap", xlab="White Obama Support", xlim=c(0,1), ylim=c(0, 0.15), 
     main="White Gender Gap Largest\nIn Contested White Counties")
symbols(agg$x, agg$y, circles=sqrt(agg$n) / 2000, 
  bg=rgb(0.5, 0.5, 0.5, alpha=0.1), fg=rgb(0.5, 0.5, 0.5, alpha=0.75), 
  add=TRUE, inch=FALSE)
L <- loess(y ~ x, weights=n, span=1, data=agg)
ord <- order(agg$x)
points(agg$x[ord], L$fit[ord], type="l")
axis(1, at=seq(0,1,by=0.25), lab=paste0(100*seq(0,1,by=0.25), "%"), lwd=par()$lwd)
axis(2, at=seq(0, 0.15, by=0.05), lab=paste0(100*seq(0, 0.15, by=0.05), "%"), lwd=par()$lwd)
abline(v=0.5)
box(lwd=par()$lwd)

dev.off()
